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Abstract 
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Lattice QCD with an even number of degenerate quark flavours is shown 
to be a limit of a local bosonic field theory. The action of the bosonic theory 
is real and bounded from below so that standard simulation algorithms can 
be expected to apply. The feasibility of such calculations is discussed, but no 
practical tests have yet been made. 
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1. Introduction 



It is well-known that the quark degrees of freedom are difficult to include 
in numerical simulations of lattice QCD. At present the probably best method 
is to integrate over the quark fields analytically and to simulate the result- 
ing effective gauge field theory using the Hybrid Monte Carlo algorithm [I]. 
This ingenious algorithm represents a significant progress compared to earlier 
techniques, but the computer resources required for QCD simulations remain 
enormous (see ref.[2], for example). 

In this paper a different computational strategy is proposed. No practical 
tests of the algorithm have yet been performed and much further work will be 
needed before the merits and limitations of the method can be assessed. The 
basic idea is to map the lattice theory to a local bosonic theory by incorpo- 
rating a suitably chosen fermion matrix inversion algorithm in the functional 
integral. Simulation methods familiar from pure gauge theories and Higgs 
models may then be applied. In particular, all degrees of freedom are treated 
stochastically and no exact inversion of the Dirac operator is required in the 
course of the simulation. 

The situation is, however, not quite so simple, because the bosonic theory 
involves a large number of complex fields coupled to the gauge field. One 
actually needs to take the number of fields to infinity if lattice QCD is to 
be reproduced exactly. To avoid excessive memory requirements, the bosonic 
theory should hence be constructed in such a way that the limit is reached 
rapidly. A concrete proposal on how to accelerate the convergence is included 
in the present paper, and a number of further technical questions concerning 
the simulation of the bosonic theory are addressed. 



2. Lattice QCD 



In the following attention is restricted to the case of two degenerate 
fiavours of Wilson quarks. Improved Wilson fermions or staggered fermions 
do not present any additional difficulties. 

We are thus concerned with gauge fields U{x,ii) and quark fields ip(x), 
tjj(x) residing on a 4-dimensional hyper-cubic lattice with spacing a and ver- 
tices X. The link variables U{x,ii) take values in SU(3) and the fermion fields 
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carry flavour, colour and Dirac indices. The lattice is assumed to be flnite 
with periodic or anti-periodic boundary conditions (boundary conditions as 
required for the QCD Schrodinger functional [6] are also tolerated). 
The action of the theory is 

S[U,^,ij] = 5jC/] + a4^V^(a;)(i? + m)V'(a;), (2.1) 

X 

where 5'p[?7] denotes the Wilson plaquette action, m the bare quark mass and 

3 

^ = i E i^A. + ^A.) - (2.2) 

At=0 

the lattice Dirac operator. Hermitean 7-matrices are assumed here and the 
covariant forward difference operator is given by 

V h'(Ij{x) = — [U{x, i^)tlj(x + afi) — '<p{x)] , (2-3) 

with fi being the unit vector in the positive direction. The backward differ- 
ence operator V* is equal to minus the adjoint of V^. 

After integrating over the quark flelds, the effective gauge held distribu- 
tion becomes 

P,ff[C/] = l[det(i? + m)f e-^»[^l. (2.4) 
The normalization constant Z (the partition function) is deflned such that 

yD[C/]Peff[C/] = 1, (2.5) 

where one integrates over all gauge flelds JJ and D[C/] denotes the usual SU(3) 
invariant product measure. Note that it is the square of the Dirac determi- 
nant which enters here, because two degenerate flavours of quarks have been 
assumed. Our aim in sect. 3 will be to flnd an exact representation of Pefr[C^] 
in terms of a local bosonic fleld theory. 

The lattice Dirac operator is a sum of an anti-hermitean and a hermitean 
part. If we deflne 75 = 70717273, it is however easy to show that 

[75(i? + m)]^ = 75(i? + m). (2.6) 
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The operator 75(1) + m) thus has a complete set of eigenvectors with real 
eigenvalues. In particular, the quark determinant det(_D + m) is real and the 
distribution (2.4) defines a probability measure. 

Without loss of generality the quark mass m may be chosen such that 
the hopping parameter k = (8 + 2am)~^ is non-negative. One may then prove 
that 

1175(1) + m)|| < M, M = 8/a+m. (2.7) 
In the following it will be convenient to work with the normalized quark matrix 

Q = -f5{D + m)/M, (2.8) 

which has eigenvalues between —1 and 1. 



3. Transformation to a local bosonic theory 

In the first step of the transformation we need to choose a polynomial 
P(s) which approximates the function 1/s in the interval < s < 1. For 
example, we may take 

n 

P(,) = ^(l_,)fc, (3.1) 

k=0 

which obviously satisfies 

lim P{s) = 1/s for < s < 1. (3.2) 

This particular polynomial is not the best choice, but to explain the transfor- 
mation it is helpful to start with a simple case. A better polynomial will be 
discussed in sect. 4. 

Since all eigenvalues of are between and 1, the matrix P(Q^) con- 
verges to the propagator In particular, 

detg2= \det P(Q^)]~\ (3.3) 

n^co 

Note, incidentally, that any polynomial approximation P(s) to the function 
1/s defines an algorithm to solve the linear system Q'^ip = 1]. One simply 
evaluates P(Q^)r] using some recursive scheme as in refs.[3-5]. 
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The roots Zk, A; = 1, . . . , ra, of the polynomial (3.1) are given by 



^fc = l-exp|2--^|. (3.4) 

For even n they come in complex conjugate pairs with non-zero imaginary 
parts. The polynomial may then be written in the manifestly positive form 

n 

nQ') = X{ [{Q-i^k? + vi], (3.5) 

k=l 

where the jikS and z/^'s are determined through 

fJ'k + ivk = V^k, Vk > 0. (3.6) 

Every value of Vk occurs twice and the corresponding jikS have equal magni- 
tude and opposite sign. 

The bosonic theory alluded to above couples the gauge field to n complex 
fields (jjkix) with colour and Dirac indices but no fiavour index (ra is taken to 
be even in the following). The action is 

n 

Sb[U,4>] = Sg[U] + ^ ^ {|(g - i^k)Mx)\' + '^l\Mx)\'} , (3.7) 

X k=l 

and from the above we now infer that 

PeAU] = lim ^ / D[<^]D[<^t]e-s.[t/,0]^ (3.3) 

where Zb denotes the partition function of the bosonic system, 

Zb = /D[C/]D[(^]D[(^t]e-s.[c/,0]. (3.9) 



It should be emphasized that the action Sb[U,(l)] is local and non-negative. In 
particular, the Gaussian integrals in eq.(3.8) are well-defined. 

The construction of the bosonic theory can be carried out in exactly the 
same way for almost any polynomial approximation P(s) to the function 1/s. 
One only requires that eq.(3.2) holds and that none of the roots Zk of the 
polynomial are real. 
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4. Chebyshev acceleration 



For the numerical simulation of the bosonic theory, n is set to some large 
but finite value, depending on the lattice size and the physical conditions. One 
must then make sure that the systematic error arising from the finiteness of 
n is negligible compared to the statistical errors. Since one cannot afford to 
take n to arbitrarily large values, it is important to achieve rapid convergence 
by optimizing the polynomial P(s). 

A similiar optimization problem occurs when solving systems of linear 
equations through an iterative procedure (see chapt. 4 of ref.[5] for a partic- 
ularly readable account). It is well-known that the convergence of some of 
the algorithms can be significantly accelerated using Chebyshev polynomials 
[3-5]. While the situation here is slightly different, the basic ideas carry over 
and lead to the optimized polynomial P(s) defined below. 

4-1 Convergence criteria 

The rate of convergence of a matrix inversion algorithm depends on the con- 
dition number p of the matrix. For a positive hermitean matrix, p is equal 
to the ratio of the largest to the smallest eigenvalue. The condition number 
of the fermion matrix is roughly proportional to the square of the inverse 
lattice spacing and is hence rather large in the cases of interest. 

The rate of convergence of the algorithm derived from the polynomial 
(3.1) can be determined as follows. First note that 



Exponential convergence is thus guaranteed for all spectral values s > 0, but 
the rate is rapidly varying and goes to zero when s becomes small. If we 
assume that the largest eigenvalue of the quark matrix is close to 1 (which 
is plausible since Q has been normalized), the condition number p will be 
approximately equal to the inverse of the smallest eigenvalue of . We then 
conclude that P(Q^) converges to exponentially with a rate close to 1/p. 

Good inversion algorithms achieve much better rates, equal to 4:/y/p or 
at least 2/y/p [5]. Our aim in the following will be to construct a polynomial 
related to such an algorithm. 

A further criterion for a good choice of polynomial is obtained when con- 
sidering the force acting on the gauge field. The force is proportional to 
the variation ^5'eff[C^] of the action of the effective gauge theory under small 




for < s < 1. 
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changes of the link variables. After replacing detQ^ by [detP(Q^)] ^, the 
variation of the action is given by 

SSMU] = SSg[U] + Tr {SQ'P'{Q')/P{Q')} , (4.2) 

where P'(s) denotes the derivative of P(s). The force, or some integrated 
form of it, drives the stochastic evolution of the gauge field in numerical simu- 
lations of QCD. We should, therefore, also make sure that the ratio P'(s)/P(s) 
converges rapidly to —1/s. 

4.2 Chebyshev polynomials 

For any real number u between and 1 an angle 9 may be defined through 

cos6' = 2m-1, < 6* < tt. (4.3) 
The (modified) Chebyshev polynomial T*{u) of degree r is then given by 

T;(u) = cos{re). (4.4) 

In particular, 

T*(u) = 1, T*(u) = 2u-1, (4.5) 
and the higher-order polynomials may be worked out recursively using 

uT:{u) = i {t:^,{u) + T:_,{u) + 2T:{u)} , r > l. (4.6) 

From the definition (4.4) it is obvious that 

sup \T:{u)\ = 1, (4.7) 

0<M<1 

and the relation 

^\^^-^^]=4T:(u) (4.8) 
du [ r + 1 r-lj '^^^ ^ ^ 

is also obtained with little effort. 

Chebyshev polynomials are uniformly bounded in the interval < m < 1. 
Outside this range they are rapidly and monotonically increasing in magnitude. 
For M < 0, for example, we have 

r;(^i) = (-l)'^cosh(rx), (4.9) 
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where x > is determined through 



coshx=l-2M. (4.10) 

4.3 Definition of P(s) 

We first seek a polynomial approximation to the function 1/s in the interval 
e < s < 1, where e is an adjustable parameter satisfying 

< e < 1. (4.11) 

To this end the fit range is mapped to the interval < m < 1 through the 
transformation 

s ^ M = (s - e)/(l - e). (4.12) 
We then consider the polynomial 

JS±lM_5LlW) (4.13, 

^ [ n+1 ra-lj ^ ^ 

and choose the constant p such that 

i2(0)=-l. (4.14) 

As will become clear below, R(s) plays the role of an error term and was chosen 

with care so as to fulfill the criteria discussed in subsect. 4.1. Explicitly one 
finds 

fcoshffra + 1)0;) coshffra — 1)0;) I 
\, , ' 4.15 
n + 1 n — 1 

where a > is determined by 

coshQ! = (l + e)/(l-e) (4.16) 

(here and below we assume that n is even). 

We now define a polynomial P(s) of degree n through 

Pis) = [l + Ris)]/s. (4.17) 
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Note that 1 + -R(s) vanishes at s = and is therefore divisible by s. From the 
definition of and the uniform boundedness of the Chebyshev polynomials 
it is clear that 

|P(s) - < -^i^L for e < s < 1. (4.18) 
Furthermore, since the constant p goes to zero exponentially at large ra, 

71 

— — e-'^«{l + 0(l/n)}, (4.19) 
smh a 

we conclude that P(s) approximates 1/s in the range e < s < 1, with an 
absolute relative error which is uniformly bounded and exponentially small. 

The polynomial in fact converges to 1/s for all s between and 1, although 
with a gradually decreasing exponential rate when < s < e. More precisely, 
the rate in this range is equal to a — x, where x is determined through eqs.(4.10) 
and (4.12). The maximal value, 

4o In 
^'»' = TT7™"l»«)= 0^7)1^^(1 + 0(1/..)}. (4.20) 

diverges linearly with ra, with a large coefficient if e is small. 
4.4 Properties of P'(s) 

As discussed in subsect. 4.1, one also requires that P'(s)/P(s) is a good ap- 
proximation to the function —1/s in the range < s < 1. 

The particular form (4.13) of the polynomial R(s) was chosen so that 

R'(-^) = (4-21) 

[cf. eq.(4.8)]. The error term in the relation 

P'(s)/P(s) = -1/s + R'(s)/[1 + R(s)] (4.22) 

(which one deduces straightforwardly from the definition of Pis)) is hence 
uniformly bounded in the range e < s < 1. In particular, at large n the bound 
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is obtained. 

When s is between and e, the ratio P'(s)/P(s) continues to converge 
to — with a reduced exponential rate equal to a — x- 

4.5 Computation of the roots Zk 

To write down the action of the bosonic theory derived in sect. 3, one needs to 
know the roots Zk of the polynomial P{s). It does not seem possible to obtain 
them in closed analytic form. They can always be computed numerically, of 
course, to any desired precision. Since the degree n can be quite large, it 
is however not advisable to use standard library routines for this task. A 
better way to proceed is described here, and we shall also obtain some useful 
qualitative information on the distribution of the roots in the complex plane. 

Since we are mainly interested in the cases where n is large, we may 
without loss assume that < < The polynomial P{s) can then be shown 
to be positive for all s and its roots Zk thus come in complex conjugate pairs 
with non-zero imaginary parts. 

A systematic expansion of the roots in powers of 1/ra may be deduced as 
follows. First note that the ra + 1 solutions of the equation 

1 + R{s) = (4.24) 

are s = and s = Zk, k = 1, . . .,ra. Now when s has a non-zero imaginary 
part, there is a unique complex number w such that 

l{w+l/w) = 2u-l, \w\>l. (4.25) 

In terms of this variable, eq.(4.24) can be rewritten in the form 

yjn+l ^ (4_26) 

where r > and f(w) are defined through 

r = {2(n+l)/pY^^''+^\ (4.27) 
f(w) = {1 - [(n + l)/(n - l)](w-^ + w-^'') + w-^'^-^y^ . (4.28) 

We are only interested in the solutions Wk of eq.(4.26) with non-zero imaginary 
part and magnitude greater than 1. Taking the n + I'th root of the equation. 
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Fig. 1. Roots of P(s) for n = 30 and e = 0.1. The thick line 
represents the range < s < 1 of the spectrum of the quark matrix . 

it is then straightforward to show that these are given by 

Wk = -Vk {1 - [{n + l)/(n - 1)]^,-'}"'^^'^+'^ {1 + 0(l/n2)} , (4.29) 
where Vk is defined through 

r 27rA; 1 , , , 

Wfc = rexp 1 2— j- , A; = l,...,ra. (4.30) 

This result translates to an asymptotic expression for the roots Zk through 

zk = i(l + e) + i(l - + (4.31) 

The approximation obtained in this way is quite accurate in most cases of 
interest. 
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For the numerical simulations one requires the roots to machine precision. 
This is now easy to achieve by applying a few iterations of Newton's algorithm 
to eq.(4.24), starting from the large n values of the roots. 

As shown by fig. 1, the roots tend to lie on an ellipse surrounding the 
spectral interval < s < 1. This pattern can be understood from the large n 
expressions derived above. When only the leading terms are kept, the formulae 
reduce to 

1 / 1 / \ f 27rA; . . . 2'iTk ^ 

zi. = Tril + s) — Tril — s) < cosh a COS h t smh a sm > . (4.32) 

2Vt;2V ;| n+1 n+lj ^ ' 

The smaller axis of the ellipse thus has a length equal to ^(1 — e) sinh a. As n 
increases the ellipse does not change very much, but the roots get denser with 
a separation between neighbours proportional to 1/ra. 

Summary 

The polynomial P(s) defined in this section depends on the degree n and a 
parameter e. For fixed e and n oo, it converges exponentially to 1/s, for all 
spectral values s in the range < s < 1. The rate of convergence is equal to 

a = 2v^+ 0(e^/2) (4.33) 

if s > e, and falls continuously to zero when s < e. 

We have also noted that the roots of P(s) come in complex conjugate pairs 
with non-zero imaginary parts if n is sufficiently large. The polynomial thus 
has all the properties required for the transformation to the bosonic theory to 
work out. It is a much better choice than the simple polynomial (3.1), because 
the rate of convergence at large n is constant over an adjustable range of s. In 
particular, by tuning e we can avoid that most of the numerical effort goes into 
decreasing the approximation errors in places where they are already small. 

We finally remark that the matrix inversion algorithm derived from this 
polynomial has an exponential rate of convergence close to 2/y^ if we set 
e = 1/p (where p is the condition number of the quark matrix). 
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5. Miscellaneous remarks 



a. Interpretation of the fields (j)k. The dominant contribution to the integral 
over (j)k in eq.(3.8) comes from a region in field space where \\{Q — i-ik)4'k\\ is 
not much greater than //fc||<?!>fc||. With high probability the field (j)k is hence 
found in the linear subspace spanned by all eigenstates of the quark matrix 
with eigenvalues in a range of width Vk (or a few times Vk) around jik- Note 
that Vk is of order ^/e and so is much smaller than 1 in the cases of interest. 
4>k may, therefore, be thought of as representing the contribution of the quark 
modes in a narrow spectral interval to the action of the effective gauge theory. 
Considering fig. 1 it is clear that the whole spectrum of the quark matrix is 
covered in this way. 

b. Choice ofn and e. As already noted before, the bosonic theory is expected to 
be accessible to standard numerical simulation techniques. In such calculations 
n should be greater or equal to some minimal value n^i^ to guarantee that the 
systematic effects stemming from the finiteness of n are negligible compared 
to the statistical errors, n^i^ depends on e, the lattice parameters, the desired 
level of precision and the quantities to be computed, e should obviously be 
tuned so as to make n^i^ is as small as possible. 

Extensive empirical tests will be needed to determine n^i^ and the opti- 
mal value of e in any given situation. We may, however, obtain some insight 
into the problem through the following qualitative argumentation. 

It is quite clear that the spectral interval e < s < 1, where uniform con- 
vergence is guaranteed, should contain most eigenvalues of the quark matrix 
. A crude free field estimate for the number of levels with s < e is 



where V denotes the lattice volume and the quark mass has been set to zero 
for simplicity. On a IG** lattice, for example, a value of e around 0.001 is 
required for at fixed physical conditions would grow to be of order 1. If 
we then take n = 380 or so, the ratio P' (s)/ P(s) approximates — 1/s with an 
absolute accuracy better than 10~^ in the range e < s < 1 [cf. eq.(4.23)]. 

These figures are perhaps a bit pessimistic and they should in any case 
not be taken too seriously. But they suggest that it may be necessary to set 
n to values as large as a few hundred. 
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c. Storage requirements. On a IG** lattice the memory size required to store 
100 fields (j)k is about 157 Mwords. Although modern parallel computers tend 
to have a core memory sufficiently large to cope with such amounts of data, 
other ways of dealing with the storage problem exist. 

On computers with fast I/O channels, for example, the field (j)k may be 
read in from an external storage device only when it is updated. This works 
out because (j)k does not couple to any field other than the gauge field (which 
is held in the main memory of the computer). Its contribution to the force 
acting on the gauge field can be computed at same time and may be added 
to the total force (which is also held in memory). After all (^^'s have been 
refreshed in this way, one may then proceed to update the gauge field. At this 
point there is no need to refer to the fields stored externally, because the force 
is already known. 

d. Critical slowing down and scaling. The most efficient simulation algorithms 
for pure gauge theories are based on the idea of over-relaxation (for a recent 
review and references see ref.[7]). It is likely that similar methods apply to 
the bosonic theory discussed here. A further acceleration of the convergence 
of the algorithm may perhaps be achieved by adopting an update cycle, where 
the gauge field and the slow modes (j)k (those with a small value of Vk) are 
visited more often than the other modes. 

If we assume that an algorithm with a dynamical critical exponent close 
to 1 can indeed be found, the computational effort at fixed physical conditions 
is expected to increase roughly like a~^, when the lattice spacing a is made 
smaller. At the same time the memory size required goes up approximately like 
a~^ . These estimates take into account that e must be scaled proportionally 
to a^ to guarantee that most eigenvalues of the quark matrix remain in the 
spectral interval e < s < 1. To preserve the level of accuracy it is then 
necessary to increase n proportionally to a~^ . 

It goes without saying that this argumentation is speculative and needs 
to be confirmed by empirical investigations. 

e. Correlation functions. From a given ensemble of gauge field configurations, 
simulating the effective gauge field distribution Peflf[?7], the meson and baryon 
correlation functions can be obtained by computing the quark propagator 
using some iterative method such as the conjugate gradient algorithm. This 
may well be the best way to proceed, but it is nevertheless interesting to note 
that the correlation functions have an exact bosonic representation. 

To deduce this representation we first observe that the correlation func- 
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tions of interest can be rewritten as the expectation value of a sum of products 
of the bilocal operator 



OaA,l3B(x,y) = ^ aA(x)llj i3Biy) 
/ = 1,2 



(5.2) 



(/ denotes the flavour index, a, (3 the colour and A, B the Dirac spinor indices). 
A generating functional for such expectation values is obtained by adding 



to the QCD action (2.1), where / is a suitable source held. The transformation 
to the bosonic theory now works out essentially as before, except that Q is 
shifted by the source term. Subsequent differentiation with respect to the 
source then yields the desired bosonic representation of the hadron correlation 
functions. 

It is worth pointing out that one ends up with correlation functions of local 
bosonic operators. The characteristic fermionic signs are taken into account 
when rewriting the product of the hadron operators in the form of sums of 
products of OaA,/3Bix,y). 



The local bosonic theory discussed in this paper provides a new starting 
point for numerical simulations of lattice QCD. Empirical studies are now 
needed to determine the computational cost of such simulations. 

As a matter of principle one might object that the method involves an 
additional source of systematic error, because the number n of bosonic flelds 
coupled to the gauge held cannot be taken to inflnity in practice. In many 
respects a flnite value of n amounts to an infrared cutoff in the quark sector 
and thus plays a role similar to the flnite extent of the lattice. To be fair one 
should also say that all known simulation algorithms for dynamical fermions 
are affected with similar systematic errors. In the case of the Hybrid Monte 
Carlo algorithm, for example, one relies on an iterative scheme to compute 
the inverse of the quark matrix. Since only a flnite number of iterations can 





x,y a, (3 A,B 



6. Conclusion 
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be performed, one then needs to show that the ensuing error has a negligible 
effect on the simulation results. 

I would like to thank Rainer Sommer, Peter Weisz and UUi Wolff for 
discussions and critical comments on a first draft of this paper. 
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